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Abstract 

In this paper we consider the single patch pseudo-spectral scheme for tensorial 
and spinorial evolution problems on the 2-sphere presented in [3,4] which is based 
on the spin-weighted spherical harmonics transform. We apply and extend this 
method to Einstein’s equations and certain classes of spherical cosmological space- 
times. More specifically, we use the hyperbolic reductions of Einstein’s equations 
obtained in the generalized wave map gauge formalism combined with Geroch’s sym¬ 
metry reduction, and focus on cosmological spacetimes with spatial S^-topologies 
and symmetry groups U(l) or U(l) x U(l). We discuss analytical and numerical 
issues related to our implementation. We test our code by reproducing the exact in¬ 
homogeneous cosmological solutions of the vacuum Einstein field equations obtained 
in [7]. 

1 Introduction 

For many interesting problems, in particular in general relativity, spherical topologies 
or play an important role. In the context of cosmological models the spherical 
Friedman-Robertson-Walker models, the Bianchi IX or the Kantowski-Sachs models are 
particularly important examples. The main difficulty for the numerical (and analytical) 
treatment of spherical manifolds is the fact that these manifolds cannot be covered 
globally by a single regular coordinate patch, and therefore the coordinate description of 
any smooth tensorial quantity inevitably breaks down somewhere. In the literature this 
problem is often referred to as the pole problem since in standard polar coordinates for 
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the 2-sphere these issues appear at the poles. Many approaches have been tried to 
deal with this issue, see for instance [22,30] and references therein. In earlier work [3,4], 
we have presented a numerical framework which applies to situations which involve the 
2-sphere. The main idea of this approach is to implement and extend the algorithm 
introduced by Huffenberger and Wandelt (HWT) in [29] to compute a transform for 
functions of given spin-weight s on the 2-sphere in terms of spin-weighted spherical 
harmonics. The concepts of the spin-weight, the so-called et/i-operators and of spin- 
weighted spherical harmonics were introduced originally in [36] and shall be reviewed 
in Section 2.3 below. As a consequence of this formalism, our code is (pseudo-)spectral 
in space; time evolutions are performed with the method of lines and standard ODE 
integrators (see below). We also point the reader to alternative implementations of this 
and similar formalisms in [2,10,25]. 

In our earlier work [3,4], we have studied simple evolution problems, like the 2 -|- 1- 
Maxwell and 2 -I- 1-Dirac equation, on fixed S^-backgrounds as test applications for our 
numerical infrastructure. The main motivation for this paper now is to apply the same 
formalism and numerical infrastructure to the much more complicated situation of the 
full Einstein equations. In this context, 2-spheres arise in a very natural way. In the 
asymptotically flat setting for example the spatial manifold M x is often considered 
which allows to address the spherical character of the far zone of the radiation field. In 
the cosmological setting, which we are interested in here, we can find S^-topologies as 
a consequence of Geroch’s symmetry reduction [24] (see Section 2.1) when the original 
spatial manifold has a symmetry. For example this was the basis for the work by Moncrief 
in [34] and for subsequent papers, and for the work in [6,7] which shall play a particularly 
important role in Section 5. Here the spacetimes of interest have spatial S^-topologies 
and the metrics have a certain spacelike symmetry such that Geroch’s reduction yields 
the spatial manifold The 3 -|- 1-vacuum Einstein equations thereby become 2 -|- 1- 
coupled Einstein-scalar field equations. All of this is explained in Section 2.1. Notice that 
Geroch’s symmetry reduction has also been used to obtain axially symmetric reductions 
of Einstein’s equations in the asymptotically flat case; see for example [11,40]. 

The extraction of suitable evolution and constraint equations from Einstein’s equa¬ 
tions is essentially the same problem, both in the original 3-1-1 and in the reduced 2 -|- 1- 
case. We use the generalized wave map formalism [21] (also called generalized harmonic 
map formalism or simply wave/harmonic map formalism) which can be understood as 
a covariant version of the more familiar generalized wave/harmonic formalism; the lat¬ 
ter was first introduced in [21] in order to generalize the original harmonic/wave gauge 
considered in [18]. We summarize the wave map formalism in Section 3.1. It turns out 
that in combination with the spin-weight formalism all singular terms caused by the 
singular polar coordinate chart of the 2-sphere can be completely eliminated. This had 
already been observed for the simpler equations considered in [3,4]. Notice that gener¬ 
alized wave gauges have been used extensively in the literature in various contexts, see 
for instance [23]. 

The numerical results in this paper are obtained using the spin-weighted spherical 
harmonics transform in [3,4]. The underlying Fourier transform is 2-dimensional as it 
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applies to functions defined on the 2-dimensional manifold Sometimes however it 
is interesting to restrict to special classes of functions on and therefore to derive 
a specialized, but more efficient version of this transform. In our context we shall be 
interested in functions on which are invariant under rotations around an axis (in 
M^), i.e., functions which do not depend on the azimuthal angle ^ in standard polar 
coordinates. For such functions, the 2-dimensional transform is inefficient. In this paper, 
we therefore also present an efficient implementation of a 1-dimensional variant of this 
transform which applies to such axially symmetric functions on The complexity 
0{L^) of the 2-dimensional transform is thereby reduced to the complexity 0{L‘^), where 
L is the band limit of the functions on in terms of the spin weight spherical harmonics. 
We call this transform the axially symmetric spin-weighted transform. It will be discussed 
in detail in Section 4. 

Finally, Section 5 is devoted to a test application of our numerical approach. We 
discuss different error sources and how they arise in our implementation. We also study 
the evolution using the areal gauge and the generalized wave map gauge. 

2 Preliminaries 

2.1 Geroch’s symmetry reduction 

In this section, we give a quick overview of Geroch’s symmetry reduction [24]. Let 
M = M X S be a globally hyperbolic 4-dimensional spacetime endowed with a metric 
Qab of signature and a global smooth time function t whose level sets are 

Cauchy surfaces homeomorphic to S. We denote the hypersurfaces given hy t = to for 
any constant to by Each is homeomorphic to E. To begin with, let be a 
smooth space-like Killing vector field on M induced by the smooth effective global action 
of the group U(l). We suppose that is everywhere tangent to the hypersurfaces Sj 
and define 

^■.= gabe^\ ^a:=eabcde^^^^ ( 2 . 1 ) 

as the norm and the twist of respectively. The operator 0 is the covariant derivative 
compatible with the metric gab- Notice that by the Frobenius theorem, the field is 
hypersurface orthogonal if and only if fla = 0. We also define 

hab Qab ^Cafb- (2-2) 

if 

It turns out that for vacuum spacetimes (M,gab) with any cosmological constant A the 
1-form rid is closed. This fact allows us to introduce a local twist potential u> so that 
In fact, tD is a global potential if M is simply connected as we shall always 

assume. 

Let S be the set of orbits of on M and consider the map 

TT : M ^ 5, 
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where tt maps every p £ M to the uniquely determined integral curve of through p. 
The requirement that tt is a smooth map induces a differentiable structure on S, and 
hence it can be considered as a smooth manifold. Since C^hab = 0, there is a unique 
smooth Lorentzian metric on S which pulls back to hab along vr. We write this metric 
on S as hab- For the same reason, there are also unique functions V', cu on S' which pull 
back to V’ and Co. It turns out that Einstein’s field equations with cosmological constant 
A for {M,gab) imply the following set of equations for {S,hab,'4’,^) where 

hab ■— '^’hab- (^•3) 


We call this system the Geroch-Einstein system (GES)^ and it reads 


VaV“u; 


h^ab 


4 (V„V’V“V’ - V„a;V“a;) 

Ip 

Ip 


Eab + 


iP 


2A, 


(2.4) 


with 

Eab = ^ bip +bOj) , (2.5) 

where Rab is the Ricci tensor associated with hab- In fact, the equations for ^p and ui 
imply that 

Tab ■= Eab - -habE, ( 2 . 6 ) 

is divergence free. It thus plays the role of the energy momentum tensor associated with 
the two scalar fields ip and uj in the 2 + 1-dimensional spacetime S with metric hab- As a 
result, the GES can be interpreted as the equations of 2-|-l-gravity coupled to two scalar 
fields Ip and oo governed by wave-map equations. 

Suppose for a moment that ip, oo and hab are known and that they satisfy Eq. (2.4). 
Then we can reconstruct a solution (M, gab) of the 3 -|- 1-dimensional Einstein vacuum 
equation with cosmological constant A as follows. It turns out that as a consequence of 
the above equations, the 2-form 


^ab — 


2 ^ 3/2 


on S is curl-free, i.e.. 


^[a*af)c] O' 


We can pull this quantity back to a 2-form ex on M which is also curl-free. This means 
that there exists, locally, a 1-form pa on M such that 


^We have used Va to denote the covariant derivative operator associated with hab- Indices are lowered 
and raised by hab- 
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notice that it is irrelevant here that the Levi-Civita covariant derivative T)a is not yet 
known at this stage due to the antisymmetrization. In fact, the 1-form rja is uniquely 
determined by this equation up to the gradient of a smooth function / and we can use 
some of the freedom in choosing / to set = 1- Then the metric 

gab = Kb + ’^fjaflb, 

where hab and are the pull-backs of /lafe/V' and V', respectively, is a solution of the 3-1-1 
Einstein vacuum equation with cosmological constant. A, and, we find 

K = 9abC’' = 'tPVa- 

So effectively, the quantities V’, bo and hab determine the metric Qab up to a total gradient 
of some function / which fixes the covariant version K of the Killing vector field 

2.2 The 2- and 3-spheres 

To begin with, we consider the manifold as the submanifold of given by x\ + x^ + 
x\ + x\ = 1. We can introduce Euler coordinates { 6 , Ai, A 2 ) on 

xi = cos - cos Ai, X 2 = cos - sm Ai, 

X 3 = sm - cos A 2 , X 4 = sm - sm A 2 , 

where 6 G ( 0 , tt) and Ai, A 2 G (0,27r). Alternatively, we use coordinates { 6 ,pi,p 2 ) (which 
are also referred to as Euler coordinates) with 9 as above and 

Ai = (pi + ^ 2 ) 72 , A 2 = (pi -/32)/2. (2.7) 

Clearly, both sets of Euler coordinates break down at 0 = 0 and vr. The vector fields dp-^ 
and dp 2 are smooth and non-vanishing vector fields on which become parallel at the 
poles 6 = 0,71. 

Similarly we define the manifold as the subset yf + ~ ^ of introduce 

standard polar coordinates 

yi = sin d cos ip, y 2 = sin t? sin p, = cos d. 

The Hopf map tt : ^ is 

(xi, X2, X3, X4) 1 -^ (yi, y 2 , ys) = (2(xiX3 -h X2X4), 2(X2X3 - X1X4), x? x^ - x| - xl) 

= (sin 9 cos p 2 , sin 9 sin p 2 , cos 9). 

This is a smooth map which has the coordinate representation 

TT : {0, Pi,P2) ^ {d, p) = {9, P2). (2.8) 
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Hence, with respect to our coordinates, the Hopf map reduces to a simple projection 
map. Now, is a principal fiber bundle over with structure group U(l) whose bundle 
map is the Hopf map. In fact, if M = M x and = 9“^ is assumed to be a Killing 
vector (as we will do) then S' = M x is the space of orbits and tt the corresponding 
map in Geroch’s symmetry reduction. 

In this paper, we shall employ this relationship between and for our studies 
of C/(l)-symmetric fields. Just as a side remark we also mention the case S = x 
which is a trivial bundle over If agrees with a vector field tangent to the S^-factor 
and we introduce appropriate coordinates, then the bundle map vr : x —)■ takes 

the same coordinate form as Eq. (2.8). In particular, Geroch’s symmetry reduction also 
yields the space of orbits 5 = M x S^. Hence, almost all techniques which we introduce 
in this paper, can also be applied to that case. 

2.3 The bundle of orthonormal frames over and spin-weighted spher¬ 
ical harmonics 

SO(3) is the bundle of oriented orthonormal frames over with structure group U{1). 
Given that SO(3) is double covered by SU(2) and that the latter is diffeomorphic to 
the Hopf map tt : —)■ can be identified with the bundle map. The theoretical details 

are discussed, for example, in [4]. Hence, when we start from the spatial manifold do 
the symmetry reduction as explained before and therefore arrive at the spatial manifold 

the manifold “reappears” in a different role, namely as the bundle of orthonormal 
frames. In practice this means the following: we let U be the dense open subset of 
obtained by removing the north and the south poles. The polar coordinates (i?, cp) cover 
U and the Euler coordinates {9,pi,p2) cover 7r~^{U). In particular, Eq. (2.8) holds. Let 
{rrf'^rrf') be the complex smooth frame on U defined by 

and by the complex conjugate m“. Any point p = {9,pi,p2) G in the bundle of 
orthonormal frames can then be identified with the basis of the tangent 

space evaluated at the point 7r(p) G The local section a : U —)-7r“^(f7) specified by 
any real function pi = pi{'{},(p) yields a different frame over U which is related to 
(m“,m“) by a pointwise rotation 

^ ( 2 . 10 ) 

at each point in f7. If / : [/ —)■ C is a component of a smooth real tensor field on 
with respect to the frame (m“,m“), the function • (/o tt) on 7r~^{U) C which 
is defined for some integer s called the spin-weight, is the corresponding component 
obtained by any frame rotation above. Any such function / is said to have the well- 
defined spin-weight s. The “standard” section and hence the “standard” frame which we 
shall use without further notice in the following is given by pi{'9,ip) = 0. We shall not 
distinguish between the original function f on U and the corresponding function / o tt 


6 



on the range of the standard section in the bundle of orthonormal frames. When we 
interpret a function / on with well-defined spin-weight s as the function • (/ o tt) 
on 7 r“^( 17 ) C we are able to replace singular frame derivatives on by regular 
derivatives along left-invariant vector fields on This yields et/i-operators 3 and 0 
given by 


S[/] := - s/cot??, (2.11) 

§[/] := du[f] ++sfcoti^, (2.12) 

for any function / on with spin-weight s. Notice that our convention differs by a 
factor \/2 from the one for and in Eq. (2.9). The function 3[/] has a well-defined 
spin-weight s -|- 1 and 0[/] has spin-weight s — 1. 

The spin-weighted spherical harmonics (SWSH) play an important role in the rep¬ 
resentation of spin-weighted functions on They form a basis of L^(SU(2)) as an 
application of the Peter-Weyl theorem to the compact group SU(2) [43]. Under certain 
assumptions, any spin-weighted function gf on can therefore be represented as an 
infinite series of SWSH 


OO I 

sf{d,(p) = EE ^Im s^lm 

1=0 m=—l 

where sEJm are the SWSH and aim the complex coefficients of the function (also called 
spectral coefficients). The standard scalar spherical harmonics are given by s = 0. By 
applying the et/i-operators to these we obtain 

{I - s){l + S + 1) s+lYimi'd, </?), 

3[sE(m(l?,V?)] = y^{l + s){l - S + 1) s-lYim{^, ^p), 

SS[sE(m(^?,</?)] = -{I - s){l + S+ 1) sYim{d,ip). 

3 Einstein’s evolution and constraint equations 

3.1 Hyperbolic reduction 

The Einstein equations (2.4) are a set of geometric partial differential equations. They 
are invariant under general coordinate transformations which implies that they are not 
automatically of any particular type when expressed in an arbitrary coordinate system. 
There exist many ways of extracting hyperbolic and elliptic subsets from these equations 
by fixing certain coordinate gauges. Here, we will use the so called wave map gauge, a 

generalization of the well-known harmonic gauge. The setup for the wave map gauge is 

discussed in detail in the appendix. 

We now introduce a general smooth frame (e“). Notice that this frame is neither 
necessarily a coordinate frame nor an orthonormal frame. The components of the 
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(3.1) 


Ricci tensor Rab with respect to this frame can be written as 

R^P = -\ + V(<,r^) + T^p{h,dh), 

where in the first term, dyh^p is the derivative of the function h^p in the direction of the 
frame vector field e“, the third term is a lengthy nonlinear expression in the components 
of hab and their first derivatives, and in the second term denotes the contracted 
connection coefficients Tp := h’-'^Tpy^j with 

vr — /IF — F*^ F 

Here, the connection coefficients of the frame are defined as (using the conventions in [42]) 


V^e“ = F^ 




e 


a 

G 


and are computed as 

flU = h^'^Tppy, ^ fiup ~ 2^ dphpy + dyhpp — dphyp + Cppy + Cypp — Cpy p ) , (3.2) 


with Cypp = hyijC"pp and 

[ep,ey\^ = C^pyCp. 

Notice that while the left side of Eq. (3.1) represents the components of a smooth tensor 
field, none of the terms on the right-hand side does this individually. In particular, the 
quantity F^ does not represent a covector field. Notice also that, in general, dpdyh^p is 
not the same as dydph^p (as it would be for a coordinate frame for which Cypp = 0). 

The non-tensorial split is not the only issue of Eq. (3.1). In addition, the second 
term destroys the hyperbolicity of its principal part. The idea is to get rid of this second 
term by defining a new tensor field 


Rab ■— Rab T (aRb) 


(3.3) 


where is the vector field defined in Eq. (A.3). The frame representation of Rab is 

Rap = -l h^'^dpdyhap + Tap{h,dh) + V,)f + V(,/^), 

where Tap is the same nonlinear expression as above. Regarded as a differential operator 
acting on hpy it has a hyperbolic principal part. 

The idea of the generalized wave map formalism is to replace the Ricci tensor Rab in 
the field equation by this new tensor Rab- We call the resulting equations the “evolution 
equations” since, under suitable conditions, these have a well-posed initial value problem 
for any choice of gauge source quantities fa and hab- The evolution equations implied 
by Eq. (2.4) are therefore 

VaV> 

Rab 
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= y (VaV’V“V'- VacaV^cu) - 2A, 

if 

2 A 

— Eab H- -hab, 

if 



with Eqs. (2.5) and (3.3). This is a coupled system of quasilinear wave equations. 

Suppose now that {hab,'ip,i^) is a solution of the evolution equations. It is a solution 
of the original equations (2.4) if = 0 (hence, if hab is in generalized wave map 
gauge) and hence Rab = Rab- Under which conditions does therefore the covector field 
2?“ vanish? The evolution equations and the contracted Bianchi identities imply the 
subsidiary system (see [39] for details) 

+ =0. (3.4) 

This is a homogeneous system of wave equations for the unknown Da- It follows that 
2 ?“ = 0 if and only if 2?“ = 0 and Va2?^ = 0 on the initial hypersurface; these conditions 
therefore constitute constraints. While the constraint 

0 = 2?'^ = /r^'^(fV-rV) + r (3.5) 

can be satisfied for any initial data hab, V' and cu by a suitable choice of the free gauge 
source quantities fa and hab, and is hence referred to as the gauge constraint, the con¬ 
straints 

= 0 (3.6) 

turn out to hold at the initial time if and only if the initial data satisfy the standard 
Hamiltonian and Momentum constraints (supposing that the gauge constraint and the 
evolution equations are satisfied). These are equations which are therefore indepen¬ 
dent of the gauge source functions. Hence Eq. (3.6) represents the actual “physical 
constraints” on the initial data for hab, V’ and oj. 


3.2 The generalized wave map gauge in the case S' = M x 

In this section, we focus on the case S := M. x and the field equations in the form 
Eq. (2.4). As before let t : S —)• M be a smooth time function on S and 

St := {t} X ~ t G M. 

We introduce coordinates {t,'d,(p) on the dense subset M x 2/ of S and define r“ = df. 
With the same choice of complex vector field on U C as in Section 2.2, we 
introduce the frame (cq, e“, e^) = (T“, m“, m“) on M x 27. The spin-weight of any function 
/ ; M X 2/ —)• C is defined in the same way as in Section 2.3, but now with respect to 
frame transformations of the form 


instead of Eq. (2.10). Therefore, the frame vector field r“ has spin-weight 0, spin- 
weight 1 and spin-weight —1. Under the above considerations, we choose the dual 

frame (cj°, w^) by 


0Ja = Vat, Lol = ^{Va'd+ i sin ttVaif) 


OJ, 


^ai 
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with spin-weight of 0, —1 and 1 respectively. The duality relation reads 

Then, the general form of a smooth metric on S is 

hab = A -h 2 (^/3 u:l^ +2S cj) ujIujI + ^ ujI^I (3.7) 

After a straightforward computation we find that almost all the quantities C^iyp intro¬ 
duced in the previous subsection are zero except 

C ‘^12 = C^ 2 i = —C '^21 = coti?. (3.8) 

The occurrence of the singular function cot is a consequence of the fact that the quan¬ 
tities C^vp are not components of a tensor field and hence do not have well-defined 
spin-weights. It is a consequence of the discussion in the previous section, however, 
that all quantities, which we eventually work with, are frame components of smooth 
tensor fields and therefore have well-dehned spin-weights without such singular terms 
— even though singular terms without well-dehned spin-weights appear in intermediate 
calculations when non-tensorial expressions are used. Since the eth-operators are essen¬ 
tially projections of covariant derivatives it is not surprising that all these terms which 
are caused by the connection coefficients related to the unit sphere will disappear when 
frame derivatives are replaced consistently by corresponding eth-operators according to 
Eq. (2.11). Indeed, we are able to demonstrate this explicitly. 

In all of what follows we choose 


hab — —^a^b + 2 ^(a^b)’ (3-9) 

as the reference metric introduced in the previous section. This is a smooth metric on 
S which represents the static cylinder with the standard spatial metric on All the 
remaining gauge freedom is then encoded in the vector field /“. We shall introduce 
quantities 

pM .= (3.10) 

f^;=r^-f^, (3.11) 

where the latter are the components of a covector field Ta which we call the smooth 
contracted connection coefficients. Thus 

_ pa _ ja 

Note that by construction, the non-tensorial quantities T^ do not contain any derivatives 
of the metric hab and we have 

Tad 

f “ftc = ( Cdbc + Cdcb - Cbcd ) 
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as a consequence of Eq. (3.2). But they contain terms proportional to coti? due to 
Eq. (3.8). All the first order derivatives of the metric hab in E^ are in Fa. We define the 
non-tensorial quantity 

t^,(/i, dh, f) := V(^f,) + T^,(/r, dh) , (3.12) 

with the same as in Eq. (3.1), and write the evolution equations as 

hP^ 

= — {dpipdai/j - dpUjd^oj) - 2A, 

hP^ 

- h^'^TP^pdpOj = — dpil^d^oj, (3.13) 

1 4A 

hP'^dpdahpu - 2T^,,(/i, dh, F) = 2V(^/^) - {dpi) d^ij) + dpUj dyUj) - -^hp„ . 

We notice that the first terms on the left hand sides constitute the principal part of this 
evolution system, i.e., quasilinear wave operators. These terms by themselves are not 
tensorial and hence give rise to singular terms (terms proportional to cotd) and terms 
which do not have well-defined spin-weights when the frame derivatives are replaced by 
eth-operators as described before. The second terms on the left hand sides cancel these 
problematic terms completely, and consequently, the left hand sides are tensorial. The 
right hand sides are tensorial already. As a result of this fully tensorial character of 
all these equations, the system of evolution equations Eq. (3.13) can now be solved by 
implementing a pseudo-spectral method based on the SWSH. 

4 Numerical implementation 

As explained earlier, we wish to implement a spectral method based on spin-weighted 
spherical harmonics to approximate spatial derivatives. A basic introduction to spectral 
methods can be found in books like [16,17,46] and references therein. For the temporal 
discretization we mainly used the Runge-Kutta-Fehlberg method except for convergence 
tests for which the explicit 4th-order Runge-Kutta method is used. We start this section 
by describing briefly the algorithm of complexity 0{L^), where L is the band limit 
of the functions on in terms of the spin-weighted spherical harmonics, to compute 
the spin-weighted spherical harmonic transforms (forward and backward) introduced by 
Huffenberger and Wandelt in [29]. Henceforth we will refer to this algorithm as HWT. 
Later, in the next subsection, we introduce an optimized version of this transform for 
the case of functions on with spin-weight s that exhibit axial symmetry (i.e., invariant 
along the coordinate vector field dp). As a result, we obtain an algorithm of complexity 
0{L?‘) which requires a low memory cost in comparison with that required by HWTs. In 
this work, we will focus on functions that exhibit axial symmetry and hence our spectral 
implementation is based on this transform. For details, improvements and applications 
of the HWTs for general functions in we refer the reader to [3,4]. We finalize this 
section by discussing our method to choose the “optimal” grid size in order to keep 
numerical errors as small as possible. 
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4.1 General description of HWTs 

To begin with, let us consider a square integrable spin-weighted function / G 
with spin-weight s. The forward and backward spin-weighted spherical harmonic trans¬ 
formations are defined, respectively, by 

saim = j fi'&, V>) sYim{"&, d^, (4-1) 

S2 

L I 

— ^ ^ ^ ^ sO^lm sYlm (??,(/?), (4.2) 

i=|s| m=—l 

where the decomposition has been truncated at the maximal mode L. Henceforth, we 
shall refer to it as the band limit. To calculate numerically the integral in Eq. (4.1) over 
a finite coordinate grid, one requires a quadrature rule and knowledge of the SWSH over 
that grid. The quadrature rule presented in [29] is based on a smooth non-invertible 
map where geometrically the poles are expanded as circles in T^. Once a quadrature 
rule on equidistant points on has been specified, we proceed to compute the SWSH 
which are written in terms of the so called Wigner d-matrices [41] by 

These matrices are easily calculated using recursion rules introduced by [45]. Adopting 
the notation := )> the Wigner d-matrices can be expressed as 

i 

iLnm = i”-" E (4.3) 

q=-l 

where n and m take integer values that run from —I to 1. Later in Section 4.2.3, we 
explain in detail how to compute the terms. The above expression allows to write 
the forward and backward spin-weighted spherieal harmonic transforms respectively as 


q=—l 

(4.4) 

f{d,ip) = 

m.n 

(4.5) 


where the matrices Imn and Jmn are computed from the standard 2-dimensional Fourier 
transforms (forward and backward, respectively) over 27r-periodic extensions of the func¬ 
tion /(■(?, (p) into T^. In general, the complexity of the outlined algorithm is O(L^). 
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4.2 Axially symmetric spin-weighted transforms 

4.2.1 The axially symmetric spin-weighted forward transform 

Let us begin by pointing out that the previous algorithm can be decomposed into two 
main tasks; namely, computation of the terms and calculation of the Imn and 

Jmn matrices by means of the 2-dimensional forward and backward Fourier transforms, 
respectively, acting on some given function /(??, ip). In what follows, we discuss in detail 
how to simplify these tasks for the case of axially symmetric functions, i.e., functions 
that only depend on the i? coordinate. 

Let us consider a square integrable axially symmetric spin-weighted function /(??) G 
with spin-weight s. Due to the (p dependence of the non-zero m modes of SWSH, 
see Eq. (4.1), we can write the function f{'d) in terms of only ^ 1 ) 0 ( 1 ?, </?). Hence, the 
forward spin-weighted spherical harmonic transform Eq. (4.1) can be written in a simple 
form as 

sa«=y /(-d) sYii^) sin'&d'ddip, 

§2 

where we have used the notation gCii = s^io and sll(i9) = Then, we rewrite 

Eq. (4.4) as 



with^ 

TT 

In := 27r j /(r?) sini? d^. (4.7) 

0 

Similarly to what is done for HWTs in [29], the number of computations required to 
obtain the spectral coefficients * 0 ; can be reduced by a factor of 2 by using symmetries 
associated with the quantities. In addition, we can introduce another reduction due 
the fact that A(jq = 0 for /-|-n = odd. This allows to reduce the number of computations 
by a further factor of 2. Therefore, we define the axially symmetric spin-weighted forward 
transform (ASET) as 


/o/ 1 ^ 

sO-i = i^Y ^ ^ho’^n^hs {n-\-= 2), (4.8) 

n=/(mod2) 

where n is a positive integer that increases in steps of two and starts at 0 or 1 depending 
on whether I is even or odd. The vector Jn is defined by 


Jn • — 


In 

In + { — l)^I(-n) 


n = 0, 
n > 0. 


^The factor 27r comes from the trivial integral over ip. 


(4.9) 
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= 0 


d — Tl 


N ^-2 

Figure 1 

The evaluation of In can be carried out by extending the function /(i?) to T = as a 27r- 
periodic function. This allows the implementation of the standard 1-dimensional Fourier 
transform in contrast to the general case of HWTs which, due to the ip dependence, 
requires a 2-dimensional Fourier transform. Now, let us define the extended function on 
T as 

^ \{-iy f{27r-^) ^>7r. ^ ^ ^ 


(4.10) 


Clearly, the vector In remains unchanged because sF{'&) agrees with f{D) on the integra¬ 
tion domain in Eq. (4.7). The function sF{'&) is chosen to be 27r periodic, hence it can be 
written as a 1-dimensional Fourier sum. However, before doing so we need to define the 
number of sampling points in T. Let us consider Fig. 1. In this diagram, the upper part 
of the circumference represents the number of samples taken for 0 < "d < vr, whereas 
the lower part shows the — 2 samples for vr < < 27r. Clearly, the subtraction by 

2 in the lower half of the circumference comes from the extraction of the poles to avoid 
oversampling. Therefore, to sample a function on T we proceed as follows. If the desired 
number of samples for a function /(i?) on is N^, then the number of samples for 
the extended function sF(t?) on T should be N!^ = — 1) and the spatial sampling 

interval will be AD = 27r/N)). Therefore, the extended function can be written as a 
I-dimensional Fourier sum by 

K/2 

sFiD) = 

k=-N^/2+l 

The substitution of this equation into Eq. (4.7) yields 


/„ = 27r ^ Fkw{k-n), 

k=-Ny2+i 


(4.11) 


where w{p) is a function Z —)■ M defined by 


2/(1—p^) p even, 

w{p) = j sinD dD = < 0 p odd, p / ±1, 

0 ±i7r/2 p = ±1. 
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By comparison with Eq. (4.11) we note that the latter is proportional to a discrete 
convolution in the spectral space. Therefore, it can be evaluated as a multiplication in 
the real space such that /„ is the 1-dimensional forward Fourier transform of 27r gF Wr 
as follows 

27r 

4 = ^ ^ exp {-inqA'd) ^F(g'A-d) Wr (g'A-i?) , 

q'=0 

where Wr{q'A-d) is the real-valued quadrature weight in T given by 

N'J2 

Wr{qA{>) = Y, w{p). 

P=-Ny2+1 

Finally, we want to emphasize that even though this way of sampling functions on T 
allows to include the value of the extended function at the poles, it yields an even 
number of modes in spectral space. Hence, we will not have the same number of pos¬ 
itive and negative modes after application of ASFT. Indeed, for the mode In ^/2 (see 
Eq. (4.9)), the vector Jy cannot be calculated since the term I_j^y 2 is not given by the 
1-dimensional forward Fourier. We avoid this issue by calculating the set of Jn terms up 
to n = N'^/2 — 1. Note that setting In '^/2 to 2 :ero does not constitute a loss of informa¬ 
tion due to the exponential decay of the spectral coefficients of the Fourier transform. 
In fact, this extra mode is in general numerically negligible and hence, it will not affect 
the accuracy of the ASFT. 

Now, in order to satisfy the Nyquist condition [17], the relation between the number 
of sampling points in T and the band limit must satisfy the inequality 

2(A^ - 1) > (2L + 1) + 1, 

where the last term on the right-hand side comes from counting the extra term without 
mirrored partner. As a result, the maximum value of the band limit for which the ASFT 
is exact is 

L = N^-2. (4.12) 

4.2.2 The axially symmetric spin-weighted backward transform 

This transform maps the spectral coefficients sO; back to the corresponding axially sym¬ 
metric function on As the inverse transform does not contain integrals, issues of 
quadrature accuracy do not arise. In a similar way as we implemented the properties of 
the 3-dimensional A(j^ term to obtain Eq. (4.8), we can write from Eq. (4.5) the axially 
symmetric spin-weighted backward transform (ASBT) as 

N'J2 

m = Y 

n=-Ar;/2+l 
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where the vector Gn is given by 


Gn := < 


L 

' E 

/=mod2 in) 


2hhl a, Al 

4 ^ sUl 


n{-s) 


nO 


if n = N^/2, 
{l+= 2). 


(4.13) 


Similar to Eq. (4.8), I increases in steps of two and starts at ^(mod2). We set Gj ^^/2 = 0 
because in the implementation of the ASET we chose In ^/2 = 0- The evaluation of 
Eq. (4.13) is carried out by a 1-dimensional inverse Fourier transform that results in 
a function sE(??) sampled on T. This function satisfies the symmetry properties in 
Eq. (4.10) where /(t?) represents the function on 0 < < tt. Thus, sT('d) 

corresponds to the extension of the function /("d) on T. 


4.2.3 Computation of the 3-dimensional 

So far the forward and backward spin-weighted spherical harmonic transforms have been 
simplified for axially symmetric functions by the implementation of a 1-dimensional 
Fourier transform instead of a 2-dimensional one as required in the algorithm HWT. 
In fact, we can simplify the computation of the terms even further. This has a 
significant effect on the efficiency of both ASET and ASBT, given that such a task takes 
around half of the execution time in practical situations. We therefore devote this section 
to discuss this issue. 

Before we explain how the terms are computed, we bring up a relevant fact for 
both ASET and ASBT. By examination of Eq. (4.8) and Eq. (4.13), we realize that we do 
not really need to calculate the complete set of A(^^ terms^ to perform the transform. 
Instead, we just need to compute up to the A(j^ term, where s is the spin-weight of 
the function that is supposed to be transformed. This yields a remarkable speed-up of 
the algorithm since in most cases s <C T. Now, based on this, we proceed to compute 
the Al^g terms implementing the recursive algorithm introduced by Trapani and Navaza 
in [45]. The recursive relations are given by the following equations 


(®) ^\o 

(b) AL 


21-1 


21 


A 


i-i 


I mi - 1 ) . z-i 

2(/-|-m)(/-|-m — 1) 


(c) 


At 


2m 


^/{l - n){l -I- n -I- 1) 


A' 


i 

(n-|-l)(m) 


' (I — n — 1){1 -|- n -|- 2) 


A 


{l-n){l + n + l) («+2)M 


where the letters “(a)”, “(6)” and “(c)” denote the sequence in which they should be 
used. We note that terms with a combination of indices outside of the correct range 
are set to 0. One way to visualize the above algorithm is by means of the pyramidal 
representation of the A(j^ terms in Fig. 2. The volume of the complete pyramid represents 

^ n and m take integer values from —I to 1. 
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the complete set of the terms. Setting the top peak of the pyramid as Aqq = 1, 
we start moving down both in the vertical direction using rule (a), and in the diagonal 
direction by (6). Thus, one can find the Aj^^ terms in the right-hand side in the front 
face of the pyramid. Then, using rule (c) repeatedly, one can find the terms behind 
the front face in order to calculate the right-hand side of the pyramid volume. If we 
need to compute the full set of AJ^^ terms, we would need to repeat this algorithm in 
order to obtain the complete right-hand side of the pyramid volume. However, we just 
need to repeat step (6) until we reach the row corresponding to I = |s| (for the given 
/-level) because we are just interested in computing the first A(^^ terms. Moreover, since 
only the Aj^^ terms with positive values of n are needed to compute both ASFT and 
ASBT (see Eqs. (4.8) and (4.13)), we apply rule (c) until we reach the column n = 0. 
The left-hand side of the pyramid volume can be obtained by applying the mirror rule 
= (—1)^“’^A^I^I (see [45]). In Fig. 3, we display a schematic representation 
of this, where the number of A(j^ terms that have to be computed are represented by 
the gray section. In this illustration we consider the collection of the Aj^^ terms for 
each l-plane. Note that the gray section is not a rectangle since we can implement the 
symmetric transposition rule A]^|^l = (—1)1^1“”'Aj^l^. In short, we will require 0{L?') 
operations to compute the A(jg terms needed for implementing both ASFT and ASBT, 
which will allow us to pre-compute the A(j^ terms for a low memory cost in comparison 
with the general algorithm for HWTs^. 

In conclusion, we have presented both the forward and backward spin-weighted spher¬ 
ical harmonic transform for the axi-symmetric case by implementing simplifications of 
the general algorithm HWTs in order to optimize them for axially symmetric functions 
in The first main simplification is the replacement of the 2-dimensional by a 1- 
dimensional Fourier transform for both the forward and backward transforms. This 
reduces the number of computations to 0(Llog2T). The second simplification lies in 
the fact that both the forward and backward transforms do not need the full set of A(^^ 

■^For L = 1024, the memory cost of AST is ^ 1MB whereas for HWT is ^ 1GB. 
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terms in the axial case. Therefore, the resulting algorithm requires 0{L?‘) operations 
for each transform. However, if we precompute the Wigner coefficients then the 
transform is only requires 0 (Llog 2 T) operations. 

These transform have been implemented in a Python 2.7 module®. Furthermore, 
the module allows to define objects that represent spin weighted functions for which an 
algebra can be defined. Hence, it can be seen not only as a set of functions, but as a 
Python environment for working with axi-symmetric SWSH. 

4.3 Choosing the optimal grid size 

Because the axially symmetric transforms are based on the Fourier transform, we ex¬ 
pect that spectral coefficients decay exponentially to zero when the band limit tends to 
infinity. Theoretically speaking, a function is described in spectral space by an infinite 
number of spectral coefficients. On the other hand, because of the machine rounding 
error®, any sufficiently smooth function is described by a finite set of spectral coefficients 
that contribute numerically to the spectral decomposition. In other words, the spectral 
coefficients with order lower than 10“^® are negligible numerically, and thus are not nec¬ 
essary for an accurate description of functions in the spectral space. Hereafter, we call 
the /-order of the last mode above order 10“^® as the optimal band limit. Consequently, 
in virtue of Eq. (4.12) the optimal band limit defines the optimal number of grid points. 
Taking a larger number of grid points than the optimal one will add unnecessary compu¬ 
tations in the transform and consequently the accuracy is reduced instead of enhanced. 
We refer to this as the sampling error. In our implementation we control this error by 
keeping the number of grid points as close as possible to the optimal case. To this end, 
we proceed as follows. Initially, we sample the initial data in a large grid. In our case 
we have chosen Nq = 1025. Then, we apply the ASFT to each function of the initial 
data and identify the highest mode which is just above the threshold 10“^®. In other 
words, we identify the optimal band limit for each function of the system. From all 
these modes we set the order of the highest mode as the optimal band limit for the 
initial data. Henceforth we will refer to this as the global optimal band limit. Using 
Eq. (4.12) we obtain the optimal number of grid points required to sample the functions 
of the system. Einally, we begin the numerical solution of the system by interpolating 
the initial data in the optimal grid. Now, we discuss how we keep the optimal grid size 
during the evolution. Eor each time step, we check the last mode of each field in order 
to observe whether they are smaller than some given tolerance. Eor this implementation 
this has been set to 10“^^. Then, if some of those modes do not satisfy the mentioned 
condition, it implies that the number of grid points is not enough for sampling some 
of the functions of the system. Therefore, we need to interpolate all the functions to a 
bigger grid. We point out that the new grid should not differ too much from the previous 

®This module can be freely downloaded under the GNU General Public License (GPL) at http: 
//gravity.otago.ac.nz/wiki/uploads/People/Axial_Spin_Weight_Functions.zip 

®In this paper, the terms “machine rounding error” and “machine precision” refer to the finite pre¬ 
cision by which numbers can be represented in a computer. We always assume that this precision is of 
the order 10“^® which corresponds to standard “double precision”. 
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one because, as we mentioned before, it could lead to too many unnecessary grid points 
and hence to larger errors. In this implementation, we decided to increase the grid by 
four points each time this is required. Using this small increment, we expect to stay 
close enough to the optimal grid and as a consequence keep a good accuracy. 

To finalize this section we point out that due to non-linearities in our evolution 
equations, some kind of filtering process is required in order to avoid the so-called aliasing 
effect. For this we use the well known 2/3-rule. For details and justification of this rule, 
see [17] and references therein. 

5 Numerical application 

5.1 Smooth Gowdy symmetry generalized Taub-NUT solutions 

It is well known that solutions of Einstein’s field equations are uniquely determined 
(up to isometries and questions of extendibility) by the Cauchy data on a Cauchy sur¬ 
face. However, there exist cases for which the uniquely determined maximal globally 
hyperbolic development [13] of the data can be extended in several inequivalent ways. 
These extensions are not globally hyperbolic and hence there are Cauchy horizons whose 
topology and smoothness may in general be complicated. Furthermore, there can exist 
closed causal curves in the extended regions which violate basic causality conditions. A 
well known example of this sort of solutions is the Taub solution [44], which is a two- 
parametric family of spatially homogeneous cosmological models with spatial topology 
Extensions through the Cauchy horizons are known as Taub-NUT solutions [37]. 

As generalizations of the Taub(-NUT) solutions, we consider now the class of smooth 
Gowdy-symmetric generalized Taub-NUT solutions introduced in [6] motivated by early 
work by Moncrief [33]. These are Gowdy-symmetric globally hyperbolic solutions of 
Einstein’s vacuum field equation with zero cosmological constant and spatial topology 
which have a past Cauchy horizon with topology ruled by closed generators. To 
cover the maximal global hyperbolic developments, the class is written in terms of the 
“areal” time function t G (0, tt) [14] and the same Euler coordinates as in Section 2.2 for 
the spatial part. In these coordinates the metrics take the form 

g = e^{—dt‘^ d0^) -|- Rq (sin^ t e^{dpi Qdp 2 )^ + sin^ 6 e^^dpl) , (5.1) 

with a positive constant Rq and smooth functions u, Q and M that depend only on t 
and 6 . A large class of such solutions of the Einstein vacuum equations were constructed 
in [6] as an application of the Fuchsian method [1]. 

5.2 A family of exact solutions 

In a subsequent paper [7], the same authors introduced a three-parametric family of 
explicit smooth Gowdy-symmetric generalized Taub-NUT solutions as an application 
of soliton methods. For this family of exact solutions, the components of the metric 
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Eq. (5.1) are given by 


where 


e^ = 




e“ = 


Rq Ue 


-M 


64c^ 1 + y ’ 


u = 4{l- x^) (1 - yf + 4c?(l + y), V = 4ci (1 - y) (1 - C3x(2 + y)), 

with X = cos "i?, y = cos t. Here ci and cs are real constants that, together with Ro, define 
particular solutions. We point out that this family of solutions contains the spatially 
homogeneous Taub solutions as the special case given by 

Cl = j + rri^ , C3 = 0, Rq = , 

with free parameters I > 0 and m G M. Inhomogeneous solutions are obtained by 
choosing any non-zero value for C 3 (see [7] for details). 

In the following we now perform the Geroch reduction described in Sections 2.1 and 
2.2 for these exact solutions. As a consequence of Gowdy symmetry, the vector field d^f, is 
a smooth Killing field of the 2-|-1 metric hah- Consequently, all the metric components of 
any Gowdy symmetric metric represented in these two coordinates are axial symmetric 
in the sense defined in Section 4 and hence the axial symmetric transform introduced 
in Section 4.2 is the natural choice for our numerical implementation discussed below. 


Before we discuss this in detail, we list the resulting formulas 

tp = i?osin^fe“, (5.2) 

dtu = -Ro ^ e^^d^Q, (5.3) 

sini7 

d^io = -Ro ^ e^^dtQ, (5.4) 

sinfl 

h = xp (e^(—df^ -|- dd^) -|- Ro sin^ 'd e^^dy^) , (5.5) 

where tp and u are the norm and twist associated with dp^ and hab the 2 -|- 1 metric. 
Next, as described in Section 3.2, we write the metric in terms of the frame (T, m,m) 
which yields 

A = sin^ t 6^’'““, (5.6) 

/3 = 0, (5.7) 

(5 = i^osin^t (e^+“ + 7?o)/2, (5.8) 

P = Rosin^t (e^+“ - 7?o)/2. (5.9) 


Henceforth, we refer to hab as the 2 -|-1 smooth Gowdy symmetry generalized Taub-NUT 
metric. 
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We notice that the quantities associated with hab are calculated from Eq. (3.11) 
by first computing the contracted Christoffel symbols of hab and then by calculating 
from Eq. (3.10) and the background metric Eq. (3.9). The results are 

fo = — cott, (5.10) 

fi = f 2 = \/2 C 3 csc^t 2 (^- 11 ) 

Here, and in all of what follows, we choose ci = 1, Rq = 2, and only vary C3. 

Eor the following it is also convenient to list the values of the metric functions at the 
time f = 7 r /2 which we shall use as the initial data for our numerical evolutions. Notice 
that we cannot use t = 0 or t = tt as initial times because the data are singular there. 
Thus, evaluating Eqs. (5.6)-(5.9) and time derivatives at t = 7 r/ 2 , we obtain’^ 


Aq — 

-4 - c| sin^ -d, 

<9tAo 

= -4c| sin^ tl, 

(5.12) 

4>o = 

— sin^ tl, 

2 

dt(po 

= c| sin^ t?, 

(5.13) 

5o = 

(2. 

4 + — sin^ -d, 

2 

dt6o 

= c| sin^ t?, 

(5.14) 

^0 = 

0, 

dtPo 

= 0. 

(5.15) 


From Eq. (5.2) and its time derivative we obtain the initial values for ipQ and dtipo 
respectively. Finally, by integrating Eq. (5.4) with respect to and setting the irrelevant 
integration constant to zero we obtain ujq. By considering Eq. (5.3) we obtain dtUjQ. The 
explicit form of these functions is 


Wo = 


V'o 


- 128(-8 + 16c3 cost?) 


(5.16) 


256 + 288C3 + 3C3 — 512c3 cos — dcg (—56 + Cg) cos 2il + Cg cos M 
8(l + icisin2il) 


n 2 ’ 


(1 — 2c3 cos 11)2 + (^ + 1^3 

128 (64 + 64c| cos^ -d — cos^ -d — 4c| sin^ 

+C 3 cos (—128 + 8c| sin^ + 9c^ sin^ t))) / B , 
dfipo = —64c3(128c 3 cos^ ?? — 32 c 3 sin^ ?? + 4 c3 sin^-i? + c| sin® 1 ? 

+16 cos t? (—12 + 5c| sin^ il) — 24c| sin^ 2il) / B , 


(5.18) 

(5.19) 


where 


B = (32 — 64 c 3 cos -i? + 64c| cos^ -d + 8 C 3 sin^ + c\ sin^ i?)" 


5.3 Numerical error sources 

The purpose of the following subsections is to describe the numerical evolution of the 
equations (3.13) for the just discussed initial data. We shall do this for two sets of gauge 

^We have used dtgo to denote the temporal partial derivative of any function g evaluated at the initial 
time. 
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source functions. Before we go into the details in Sections 5.4 and 5.5, however, let us 
discuss possible numerical error sources which we shall refer to in our discussion of our 
numerical results below. 

Clearly the time and spatial discretization gives rise to numerical errors. In general 
it is expected that time discretization errors are larger than spatial ones thanks to the 
rapid (exponential) convergence of the latter. In order to investigate the presumably 
more significant time discretization errors we shall use two different time discretization 
schemes, the (non-adaptive) 4th-order Runge-Kutta scheme and the (adaptive) Runge- 
Kutta-Fehlberg (RKF) scheme. See [16] for details about adaptive Runge-Kutta meth¬ 
ods. Spatial discretizations shall always be based on our adaptive framework discussed 
in Section 4.3. For runs using the adaptive RKF scheme, we can therefore expect that 
all discretization errors can be made sufficiently small by choosing suitable tolerance 
parameters. 

In our numerical experiments we identify further error sources which turn out to be 
particularly severe. Recall from Section 4.3 that we choose the same band limit for all 
unknowns. However, in most of our practical examples, only few of the unknowns actu¬ 
ally require high spatial resolutions. As a consequence, many unknowns are oversampled, 
which is not only inefficient numerically, but also generates undesired numerical noise. 
The origin of this noise is that the “unnecessary” modes associated with too large band 
limits are in general not zero numerically. In fact, while they are typically of the order 
of the machine precision initially, they may grow during the evolution in particular due 
nonlinear coupling of modes. Typically, the larger the difference between the optimal 
band limit for any particular unknown and the global band limit is, the larger is this 
noise. This error is difficult to control in practice and it is quite common that once this 
noise has started to grow during the evolution it continues to grow increasingly fast. We 
measure this error by looking at the evolution of the highest modes of certain represen¬ 
tative unknowns during the evolution. The only conceivable cure of this problem would 
be to work with higher machine precisions, which would, however, significantly slow 
down the numerical runs. Our numerical infrastructure is completely based on “dou¬ 
ble precision”. We have not attempted to work with higher machine precisions such as 
“quad precision” yet. Further comments on this in the context of a different numerical 
infrastructure can be found, for example, in [2]. 

Another severe, but not fully independent numerical error is associated with the 
violation of the constraints. Recall that due to Eq. (3.4), the condition T>^ = 0 is iden¬ 
tically satisfied during the evolution if (i) the evolution equations hold exactly and (ii) 
the constraints Eqs. (3.5) and (3.6) are satisfied initially. Eor our numerical calculations, 
however, both of these conditions are violated. Let us for the sake of this argument imag¬ 
ine that the constraints are violated at the initial time, but that the evolution equations 
hold identically (i.e., we pretend that the numerical evolutions are done with an infinite 
resolution in space and time and with infinite machine precision). Then Eq. (3.4) de¬ 
scribes the (exact) evolution of the in general non-zero constraint violation quantities 
T>fj,. Since the initial data for these quantities are now assumed to be non-zero, their 
evolution is in general also non-zero. Depending on the particular properties of the evo- 
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lution system and hence of Eq. (3.4), these quantities may in fact grow rapidly during 
the evolution. If this is the case, the constraint violation error can become large very 
quickly even if it is small at the initial time, and the resulting numerical solutions of Ein¬ 
stein’s equations therefore become useless quickly. This situation cannot be improved by 
increasing the numerical resolution. In fact, this error is a consequence of the structure 
of the continuum evolution equations. Various ways to reconcile this problem have been 
proposed in the literature. One of the most promising ideas [5, 9, 26] is to introduce 
constraint damping terms, i.e., to add terms to the evolution equations (i) which are 
proportional to the constraint violation quantities (hence the solutions of the evolution 
equations for the actual case of interest = 0 are unchanged) and (ii) which, however, 
turn the surface = 0 into a future attractor for Eq. (3.4). This technique has proved 
to be quite useful to produce stable calculations for asymptotically flat spacetimes (see 
for instance [8,28,31,38]). The analytic derivation of suitable constraint damping terms 
is in general difficult and is usually done based on approximations which may only hold 
in certain regimes of the evolution (see e.g., [20]). In this paper here we work without 
constraint damping terms. Nevertheless, we remark that thanks to the close relationship 
of our formulation of Einstein’s equations with the ones used in the above references, 
similar choices of constraint damping terms are expected to be useful in reducing the 
constraint violation errors in our applications. Indeed we have already gathered some 
promising experience with constraint damping terms of the type used in [38] which we 
shall report on in a future article. 

5.4 Numerical evolutions in areal gauge 

In this section now, we shall hx the gauge freedom for the evolution equations by iden¬ 
tifying the gauge source functions with the contracted Christoffel symbols given by 
Eqs. (5.10) - (5.11); we recall that we have implicitly always assumed Eq. (3.9) as the 
reference metric and continue to do so. As common in the literature, we refer to this 
coordinate gauge as areal gauge. We shall then evolve the evolution equations (3.13) 
for the initial data given by Eqs. (5.12) - (5.19) at t = Txj2 using these gauge source 
functions. The resulting numerical solutions are given in the same coordinates as the 
exact solution and direct comparisons between the exact and the numerical solutions 
can be performed conveniently by considering the error quantity 

^(i) := maf II IIl2(S2) > 


where {t, i?) represents a frame component of the exact metric for any given t, whereas 
{t, 1 ?) represents the numerical value. The norm ||•|| 2 , 2 (§ 2 ) is approximated numerically 
by the discrete f^-norm of the grid function vector. Notice that the same spacetimes in 
the same coordinates have been constructed numerically with different methods in [27]. 
However, in contrast to our discussion here, some of Einstein’s equations turn out to be 
formally singular in the “interior” of the Gowdy square in the formulation used there 
and hence are ignored to avoid serious numerical problems. 
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Figure 4: Convergence test, C 3 = 0.2. 




Figure 5: Norm of the Killing vector, C 3 = 0.2. 


Figure 6 : Twist of the Killing vector, C 3 = 0.2. 


As a first test for our numerical implementation we present a convergence test in 
Fig. 4 for C 3 = 0.2. The evolution is carried out with the (non-adaptive) 4th-order 
Runge-Kutta scheme. The figure shows the expected convergence rate demonstrating 
that the time discretization error is dominant here. This is not surprising since at each 
t all the metric components are very smooth functions that can be resolved on the grid 
with high accuracy so long as t does not get too close to t = tt. The oversampling 
and constraint violation errors discussed in the previous subsection are small during this 
early phase of the evolution. 

Next, we replace the non-adaptive 4th-order Runge-Kutta scheme by the adaptive 
RKF method. In Fig. 5 and Fig. 6 we show the numerical evolutions of the geometric 
quantities and co for C 3 = 0.2. The numerical error in these calculations are shown in 
Fig. 7 for different values of C3. The error tolerance Tol of the RKF method is chosen to 
be 10“®. This figure suggests that the numerical errors here remain bounded for a long 
time. The larger C 3 is, however, and hence the more inhomogeneous the solution is, the 
more rapidly the numerical errors grow close to t = tt as expected. Fig. 8 indicates that 
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Figure 7: Error propagation for various values 
of cz and Tol = 10“®. 



Figure 8 : Error propagation for various values 
of Tol and C 3 = 0.3. 



Figure 9: Constraints propagation for various values of Tol and C 3 = 0.3. 


the behavior close to t = vr cannot be improved by decreasing the value of Tol. This 
suggests that close to t = tt the numerical errors are not dominated anymore by the time 
discretization error, but that one of the other error sources discussed in Section 5.3 takes 
over. Our experience suggests that in fact both the oversampling error and the constraint 
violation error are significant at late times, in particular, for larger values of C3. As a 
consequence of Eqs. (5.12) - (5.19), the required band limits for the metric components 
and their time derivatives are small, but the required band limits to resolve tpo, loq and 
their time derivatives are relatively large. This discrepancy, which we associate with the 
oversampling error, is in fact larger the larger C 3 is. As already mentioned before, the 
noise generated by oversampling indeed grows during the evolution. 

In order to measure the constraint violation error, we define the quantity 

V{t) := max || ||l 2 (§ 2 ). 

In Fig. 9, we show the evolution of this quantity for C 3 = 0.3. At late times, the curves 
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look very similar to the ones of E{t) in Fig. 8. This suggests that the constraint violation 
error contributes significantly to the total numerical error. 


5.5 Numerical evolutions in wave map gauge 

In this section we describe numerical computations for the same spacetimes as before, 
but using a different coordinate gauge. To this end, we want to choose the same initial 
data as before, but work with different gauge source functions. Both the gauge constraint 
Eq. (3.5) and the physical constraints Eq. (3.6) clearly have to be satisfied at the initial 
time. Since we do not want to resolve these complicated nonlinear PDEs, our strategy 
is to use exactly the same initial data for the values of the metric components and their 
first time derivative values, and also exactly the same initial values of the gauge source 
functions as before. In order to implement a different coordinate gauge, we then apply 
the following “gauge driver condition” during the evolution whose purpose is to rapidly 
drive the gauge source functions from their initial values fixed by the gauge constraint 
towards the target gauge source functions /^: 

U = + (5-20) 


Here the parameter q controls how rapidly the gauge is driven towards the target. The 
quantities T^lto are calculated from the initial data and are understood as functions of 
the spatial coordinates only. Notice that different gauge drivers for the generalized wave 
representation of Einstein’s equations were considered in [32]. Eq. (5.20) guarantees that 
the gauge constraint is satisfied at the initial time. As discussed at the end of Section 3.1, 
the physical constraints, even though they pose highly non-trivial restrictions on the 
choice of the initial data because they are essentially linear combinations of the well- 
known Hamiltonian and momentum constraints, turn out to not be restrictions on the 
gauge source functions. Hence it is not necessary to introduce terms in Eq. (5.20) which 
account for the first time derivative of T^ at t = to- 

We apply this idea to calculate the same spacetimes as before, but now we choose 
the wave gauge as the target gauge, which is defined by the condition = 0. Eor our 
numerical tests we choose g = 10 in Eq. (5.20). Before we present our numerical results 
we notice that it is straightforward to derive the formula 


TT 1 

Hw) = 7T + 7T log 


1 — COS t 
1 -|- cos t 


(5.21) 


which for our spacetimes relates the time coordinate t in areal gauge (used in Section 5.4) 
and the time coordinate in wave map gauge. This formula holds identically even 
though Eq. (5.20) is strictly speaking not the exact wave map gauge. However as a 
consequence of ro|to=,r /2 = 0 which follows from Eq. (5.10), the target gauge source 
function /o = 0 agrees identically with /o = 0. Eq. (5.21) is then obtained by solving 
the exactly homogeneous wave map equation for the wave map time coordinate function 
with appropriate initial conditions. Eq. (5.21) allows us to make direct comparisons 
between our results here and the results in the previous section. In particular, it reveals 
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that the wave time slices = const are the same as the areal time slices t = const 
(for different constants), and, the “singularities” at t = 0 , 7 r are shifted to infinity, in 
particular, —)■ oo for t —)■ tt. We point out however that it is not possible to derive 

a formula which relates the spatial coordinates in both gauges. This is true even if 
q in Eq. (5.20) was so large that we could consider our gauge as the exact wave map 
gauge. This is a consequence of the fact that the homogeneity of the wave equations for 
the spatial wave map coordinates is destroyed by terms given by the reference metric 
Eq. (3.9). In fact we shall demonstrate below that the spatial coordinates on each time 
slice are different in areal and wave map coordinates. 

In order to obtain a more geometric and detailed comparison of the two gauges, we 
consider the Eikonal equation following [19] 

VarV“r = -1. (5.22) 

Let r be a smooth solution of the initial value problem of the Eikonal equation with 
smooth initial data tq : Sq —)• M prescribed freely on any smooth Cauchy surface Sq 
in any smooth globally hyperbolic spacetime. The method of characteristics applied to 
this PDE allows to prove that such a solution indeed always exists at least sufficiently 
close to the initial hypersurface Eq. Eor definiteness now we restrict to the case of zero 
initial data tq = 0 for all of what follows. Fix any point p in the timelike future of 
Eq in the spacetime and consider any timelike geodesic through p (with unit tangent 
vector). Any such geodesic must intersect Eq at some point xq in the past of p. There is 
precisely one such timelike geodesic through p with unit tangent vector which intersects 
Eq perpendicularly in xq and hence the point xq is uniquely determined. The value 
t { p ) of the solution r of the Eikonal equation with zero initial data then represents 
the proper time along this timelike geodesic from xq to p. The quantity r is therefore 
a meaningful geometric scalar quantity which can be used to compare our numerical 
spacetimes, in particular when the same spacetime is calculated in different coordinate 
gauges. We proceed as follows. For initial data parameters Rq = 2, ci = 1 and C 3 = 0.1 
(see Section 5.2); 

1. We calculate the corresponding solution of Einstein’s evolution equations in areal 
gauge (in the same way as in Section 5.4) and of the Eikonal equation Eq. (5.22) 
(with zero initial data) up to t = 3. The value of the resulting r function on the 
t = 3-surface expressed with respect to spatial areal coordinates yields the dashed 
curve in Fig. 10. 

2. Eq. (5.21) implies that t = 3 corresponds to t^^) ~ 4.217. For the same initial 
data parameters as in the first step, we calculate the corresponding solution of 
Einstein’s evolution equations in wave gauge numerically (using the gauge driver 
condition Eq. (5.20) with q = 10) and of the Eikonal equation Eq. (5.22) (with 
zero initial data) up to ss 4.217. The value of the resulting r function on the 

Ri 4.217-surface expressed with respect to spatial wave map coordinates yields 
the continuous curve in Fig. 10. 
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Figure 10: Proper time comparison. Figure 11: Difference of the proper times. 



Figure 12: Comparison of the constraint violations. 


Since the t = 3-surface and the ~ 4.217-surface represent the same geometric surface 
in our spacetime and since r is a geometric scalar quantity, the value of the solution of 
the Eikonal equation on this surface should be the same function in both steps above. 
However, since this function is expressed in terms of different spatial coordinates, namely 
areal coordinates in the first step and wave map coordinates in the second step, the 
two curves in Fig. 10 are slightly different. Hence Fig. 10 can be understood as a 
representation of the difference of these two sets of spatial coordinates. This difference is 
emphasized in Fig. 11 where the two curves in Fig. 10 are subtracted directly. Intuitively, 
these two sets of spatial coordinates should agree at geometrically distinct points, namely, 
at the poles and also at the equator as a consequence of a reflection symmetry which is 
inherent to our particular class of exact solutions. Indeed the difference curve in Fig. 11 
is zero at the poles 0 = 0, tt and the equator 6 = 7r/2. 

Next, we present plots of the constraint violations in both gauges; see Fig. 12. The 
dashed curve has been calculated in areal gauge (in the first step above). The continuous 
curve has been calculated in wave map gauge (in the second step above), but has then 


28 






























































































been expressed in terms of the areal time function by means of Eq. (5.21). It is interesting 
to notice that the constraint violations are significantly smaller in wave map gauge than 
they are in areal gauge towards the end of the numerical evolution. 

Finally we comment on the fact that in wave map gauge the shift quantity /3 in 
Eq. (3.7) is a non-trivial non-zero function in contrast to areal gauge; see Eq. (5.7). When 
/3 cannot be assumed to be zero identically, the algebraic complexity of the evolution 
equations is increased dramatically. It is surprising that irrespective of this it appears 
that we get better numerical results in wave map than in areal gauge. 

6 Discussion 

The purpose of our work here was to introduce a numerical approach to solve the Cauchy 
problem for spacetimes which involve the manifold We employ a fully regular rep¬ 
resentation of the Einstein equations based on the spin-weighted spherical transform 
and the generalized wave map formalism. This allows us to account for all singular 
terms explicitly which usually arise as a consequence of the coordinate singularities of 
polar coordinates at the poles of the 2-sphere. Our numerical infrastructure is based on 
the spin-weight formalism and corresponding transforms introduced in [3,4]. We have 
extended this infrastructure so that it now provides an efficient treatment of axially sym¬ 
metric functions on the 2-sphere, reducing the complexity 0{L^) of the full transform 
to the complexity 0{L?‘). We therefore expect this method to be useful also for other 
applications in future work. We have also demonstrated the consistency and feasibility 
of our approach by means of numerical studies of certain inhomogeneous cosmological 
solutions of the Einstein’s equations. 

As another application of this method we are currently studying the critical behavior 
of perturbations of the Nariai spacetime [3,4]. In particular it is suggested that larger 
amplitudes of the perturbations, which had not been studied before, could lead to the 
formation of cosmological black holes. It would be of great interest to explore the thresh¬ 
old solutions and the expected cosmological black hole solutions as well as consequences 
for the longstanding cosmic no-hair and cosmic censorship conjectures. Other conceiv¬ 
able interesting applications where our numerical infrastructure can be applied directly 
are Robinson-Trautman solutions [15] and Ricci flow [35]. 

A The generalized wave map formalism 

Whether we want to solve the Cauchy problem for the 3 -|- 1-Einstein vacuum equations 
Gah + ^Qab = 0 or for the 2 -|- 1 equations (2.4), the first task is always to extract hy¬ 
perbolic evolution equations and constraint equations with well-understood propagation 
properties from the equation for the Ricci tensor of the unknown metric. We shall now 
briefly discuss the “generalized wave map formalism”. In most of the literature, the 
related (but not covariant) generalized wave/harmonic formalism is used. While this 
is sufficient for many applications, it is a drawback for us. In fact, for applications 
with spatial S^-topologies covered by a single singular polar coordinate system, it is far 
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more convenient to work with actual covariant quantities (i.e., smooth tensor fields). 
The reason is that frame components of smooth tensor fields on have well-defined 
spin-weights (despite of the fact that the frame itself is singular at the poles) so that 
they are expandable in spin-weighted spherical harmonics, which are globally defined 
regular ‘functions’ on the 2-sphere (even though their coordinate representation may be 
singular). It turns out that expressing everything with respect to these bases, renders 
the equations manifestly regular. 

To this end, we discuss the geometric formulation of the wave map gauge [21]. We 
consider a map <I> : M —)• M between two general smooth 4-dimensional manifolds M 
and M (or open subsets thereof) equipped with Lorentzian metrics® hab and hah- The 
map is called a wave map if it extremizes the functional 


Jm 


In coordinate charts (x^) on M and (?/“) on M we obtain the Euler-Lagrange equations 
for the coordinate representation y" = y“(x^) of 4* 




dyh Qyl 

dx>^ dx'^ 


= 0 . 


(A.l) 


Here, are the Christoffel symbols for the metric hab ia the coordinate basis on 

M, and \~\h is the wave operator for scalar functions defined by hab- This equation is 
called the wave-map equation. More details can be found in [12]. If the manifolds were 
Riemannian then the analogous equation would characterize a harmonic map between M 
and M. Let us point out that the left hand side of the equation defines a geometric object, 
namely a section in the pull-back bundle <I>*TM. This is not immediately obvious due 
to the appearance of the Christoffel symbols in the second term. However, the tensorial 
character of that term under change of coordinates in M is compensated for by the first 
term which, by itself, is also non-tensorial under such coordinate transformations. 

The generalized wave-map equation is the equation (A.l) with a non-vanishing, ar¬ 
bitrary right hand side, a section in ^*TM with coordinate representation /“ 

= -/“. (A.2) 

" dxt^ dx'^ •’ ^ ’ 

The minus sign on the right-hand side is a matter of conventions. Suppose now that 
M = M and = idM- Then (x^) and (y“) are two coordinate charts for M and (A.2) 
can be read as an equation determining the coordinate system (y“) for M by imposing 
a geometrical gauge condition. This equation is a semi-linear wave equation of M which 
has solutions near any Cauchy surface so that such a coordinate gauge always exists 
locally. 

®A11 of the following arguments also hold if M and M are n-dimensional manifolds for some arbitrary 
positive integer n and if hab is a general smooth pseudo-Riemannian (not necessary Lorentzian) metric. 
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Choosing the coordinates according to this gauge, i.e., putting and express¬ 

ing the wave operator in these coordinates yields the equation 

+ f = -r, 

where are the Christoffel symbols of the metric h on M. In this equation the 

tensorial character becomes manifest since the left hand side involves the difference of 
two connection coefficients so it gives the components of a vector field in the coordinate 
basis of the (x^). Therefore, this equation holds in any basis on M as long as we 
interpret the Christoffel symbols as the connection coefficients with respect to the chosen 
basis. Note also that this implies that imposing the equation (A.2) does not constitute 
a condition on the coordinate system (x“) but a condition on the metric components 
in their dependence on the coordinates. We define the vector field 1?“ in terms of its 
components 

:= + r- (A.3) 

So, = 0 when (A.2) is imposed. A metric hab which is restricted by 2?“ = 0 is said 
to be in wave map gauge (with respect to hab)] in Eq- (3.9) we fix a particular metric 
hab- We point out that the wave map gauge reduces to the widely used generalized 
wave/harmonic gauge characterized by = —/^ on space-times with topology 
when the Minkowski metric in Cartesian coordinates x^ is used as a reference metric hab- 
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